rm(list=ls(all=TRUE))
#set working directory

#install.packages("car")
library(car)
#install.packages("MASS")
library(MASS)
#install.packages("foreign")
library(foreign)
#install.packages("readstata13")
library(readstata13)
#install.packages("list")
library(list)

mydata <- read.dta13("Zhubajie.dta",encoding="UTF-8")

# drop obs that failed attentiveness check
mydata= subset(mydata,S2==1)

mydata$age2<-Recode(mydata$age,"18:25=1;26:35=2;36:100=3;else=NA")

mydata$Edu<-Recode(mydata$educ,"c(1,2)=1;3=2;else=NA")

mydata$pubemp<-Recode(mydata$empstatus,"c(1,2,3,6)=1;c(4,5,7,8,9,10,11)=0;else=2")
mydata$IndInc<-Recode(mydata$personalinc,"c(1,2,3)=1;c(4,5,6,7,8)=2;c(9,10,11)=3;12:19=4;else=NA")

### Create Age Dummies
mydata$agedum1<-Recode(mydata$age2,"1=1;2:3=0;NA=NA")
mydata$agedum2<-Recode(mydata$age2,"2=1;c(1,3)=0;NA=NA")
mydata$agedum3<-Recode(mydata$age2,"3=1;1:2=0;NA=NA")

### Educ Dummies
mydata$educ1<-Recode(mydata$Edu,"1=1;2=0;NA=NA")
mydata$educ2<-Recode(mydata$Edu,"2=1;1=0;NA=NA")

### Pubemp Dummies 0=private, 1=public, 3=out of labor force
mydata$public<-Recode(mydata$pubemp,"1=1;c(0,2)=0;NA=NA")
mydata$private<-Recode(mydata$pubemp,"0=1;c(1,2)=0;NA=NA")
mydata$outlabforce<-Recode(mydata$pubem,"2=1;0:1=0;NA=NA")

### Individual Income Dummies
mydata$IndInc1<-Recode(mydata$IndInc,"1=1;c(2,3,4)=0;NA=NA")
mydata$IndInc2<-Recode(mydata$IndInc,"2=1;c(1,3,4)=0;NA=NA")
mydata$IndInc3<-Recode(mydata$IndInc,"3=1;c(1,2,4)=0;NA=NA")
mydata$IndInc4<-Recode(mydata$IndInc,"4=1;c(1,2,3)=0;NA=NA")

####################### List 1: Distrust Central Leaders #######################
mydata$TrustLeaders<-mydata$Q147
mydata$TrustLeaders[mydata$ListRandom1==2]<-mydata$Q149[mydata$ListRandom1==2]

mydata$ListTreat1<-ifelse(mydata$ListRandom1==1,1,0)
mydata$TrustLeaders.d<-recode(mydata$Q218,"1=1;0=0;else=NA")

data.list.1<-na.omit(data.frame(with(mydata,cbind(TrustLeaders,ListTreat1,
    agedum1,agedum2,female,educ1,IndInc1,IndInc2,IndInc3,public,private,ccp,
    urban,managers))))

data.sens.1<-subset(na.omit(data.frame(with(mydata,cbind(TrustLeaders.d,ListTreat1,
    agedum1,agedum2,female,educ1,IndInc1,IndInc2,IndInc3,public,private,ccp,
    urban,managers)))),ListTreat1==0)

###### Test of Design Effect
ListData1.test<-na.omit(data.frame(mydata[,c("TrustLeaders","ListTreat1")]))

de.test1<-ict.test(ListData1.test$TrustLeaders,ListData1.test$ListTreat1, 
                   J=3,alpha=0.05, gms=T, pi.table=T)
print(de.test1)

################
fit.list.1 <- ictreg(TrustLeaders ~ agedum1 + agedum2 + female + educ1 + IndInc1 
            + IndInc2 + IndInc3 + public + private + ccp + urban + managers , 
            J = 3, data = data.list.1, treat="ListTreat1", method = "ml", fit.start = "lm")

fit.list.1a <- ictreg(TrustLeaders ~ 1, 
            J = 3, data = data.list.1, treat="ListTreat1", method = "lm")

fit.sens.1 <- glm(TrustLeaders.d ~ agedum1 + agedum2 + female + educ1 + IndInc1 
            + IndInc2 + IndInc3 + public + private + ccp + urban + managers, 
                  data = data.sens.1, family = binomial("logit"), x=T)

#################List 2: Critisim of Govt ######################################
mydata$CriticizeGovt<-mydata$Q151
mydata$CriticizeGovt[mydata$ListRandom2==2]<-mydata$Q153[mydata$ListRandom2==2]
mydata$ListTreat2<-ifelse(mydata$ListRandom2==1,1,0)
mydata$CriticizeGovt.d<-recode(mydata$Q220,"1=1;2=0;else=NA")

data.list.2<-na.omit(data.frame(with(mydata,cbind(CriticizeGovt,ListTreat2,
         agedum1,agedum2,female,educ1,IndInc1,IndInc2,IndInc3,public,private,ccp,
         urban,managers))))
data.sens.2<-subset(na.omit(data.frame(with(mydata,cbind(CriticizeGovt.d,ListTreat2,
          agedum1,agedum2,female,educ1,IndInc1,IndInc2,IndInc3,public,private,ccp,
          urban,managers)))),ListTreat2==0)

###### Test of Design Effect
ListData2.test<-na.omit(data.frame(mydata[,c("CriticizeGovt","ListTreat2")]))

de.test2<-ict.test(ListData2.test$CriticizeGovt,ListData2.test$ListTreat2, 
                   J=3,alpha=0.05, gms=T, pi.table=T)
print(de.test2)

################
fit.list.2 <- ictreg(CriticizeGovt ~ agedum1 + agedum2 + female + educ1 + IndInc1 + 
              IndInc2 + +IndInc3 + public + private + ccp + urban + managers, 
              J = 3, data = data.list.2,treat="ListTreat2", method = "ml", fit.start = "lm")

fit.list.2a <- ictreg(CriticizeGovt ~ 1, 
                     J = 3, data = data.list.2,treat="ListTreat2", method = "lm")

fit.sens.2 <- glm(CriticizeGovt.d ~ agedum1 + agedum2 + female + educ1 + IndInc1 
              + IndInc2 + +IndInc3 + public + private + ccp + urban + managers, 
                  data = data.sens.2, family = binomial("logit"))

########################List 3: Witnessed Corruption ###########################
mydata$CorruptionExp<-mydata$Q155
mydata$CorruptionExp[mydata$ListRandom3==2]<-mydata$Q157[mydata$ListRandom3==2]
mydata$ListTreat3<-ifelse(mydata$ListRandom3==1,1,0)
mydata$CorruptionExp.d<-recode(mydata$Q222,"1=1;2=0;else=NA")

data.list.3<-na.omit(data.frame(with(mydata,cbind(CorruptionExp,ListTreat3,
      agedum1,agedum2,female,educ1,IndInc1,IndInc2,IndInc3,public,private,ccp,
      urban,managers))))
data.sens.3<-subset(na.omit(data.frame(with(mydata,cbind(CorruptionExp.d,ListTreat3,
      agedum1,agedum2,female,educ1,IndInc1,IndInc2,IndInc3,public,private,ccp,
      urban,managers)))),ListTreat3==0)

###### Test of Design Effect
ListData3.test<-na.omit(data.frame(mydata[,c("CorruptionExp","ListTreat3")]))

de.test3<-ict.test(ListData3.test$CorruptionExp,ListData3.test$ListTreat3, 
                   J=4,alpha=0.05, gms=T, pi.table=T)
print(de.test3)

################
fit.list.3 <- ictreg(CorruptionExp ~ agedum1 + agedum2 + female + educ1 + IndInc1 
              + IndInc2 + +IndInc3 + public + private + ccp + urban + managers,
              J = 4, data = data.list.3,treat="ListTreat3", method = "ml", fit.start = "lm")

################
fit.list.3a <- ictreg(CorruptionExp ~ 1,
                     J = 4, data = data.list.3,treat="ListTreat3", method = "lm")

fit.sens.3 <- glm(CorruptionExp.d ~ agedum1 + agedum2 + female + educ1 + IndInc1 
              + IndInc2 + +IndInc3 + public + private + ccp + urban + managers, 
                  data = data.sens.3, family = binomial("logit"))

######################### List 4: Bribed Cadres ################################
mydata$Bribery<-mydata$Q159
mydata$Bribery[mydata$ListRandom4==2]<-mydata$Q161[mydata$ListRandom4==2]

mydata$ListTreat4<-ifelse(mydata$ListRandom4==1,1,0)
mydata$Bribery.d<-recode(mydata$Q224,"1=1;2=0;else=NA")

data.list.4<-na.omit(data.frame(with(mydata,cbind(Bribery,ListTreat4,
      agedum1,agedum2,female,educ1,IndInc1,IndInc2,IndInc3,public,private,ccp,
      urban,managers))))

data.sens.4<-subset(na.omit(data.frame(with(mydata,cbind(Bribery.d,ListTreat4,
      agedum1,agedum2,female,educ1,IndInc1,IndInc2,IndInc3,public,private,ccp,
      urban,managers)))),ListTreat4==0)

###### Test of Design Effect
ListData4.test<-na.omit(data.frame(mydata[,c("Bribery","ListTreat4")]))
de.test4<-ict.test(ListData4.test$Bribery,ListData4.test$ListTreat4, 
                   J=3, gms=T, pi.table=T)
print(de.test4)

################
fit.list.4 <- ictreg(Bribery ~ agedum1 + agedum2 + female + educ1 + IndInc1 + 
              IndInc2 + IndInc3 + public + private + ccp + urban + managers, 
              J = 3, data = data.list.4, treat="ListTreat4", method = "ml", fit.start = "lm")

fit.list.4a <- ictreg(Bribery ~ 1, 
                     J = 3, data = data.list.4, treat="ListTreat4", method = "lm")

################
fit.sens.4 <- glm(Bribery.d ~ agedum1 + agedum2 + female + educ1 + IndInc1 + 
              IndInc2 + +IndInc3 + public + private + ccp + urban + managers, 
                  data = data.sens.4, family = binomial("logit"))

################################################################################
###################### Function: Predicted Probabilities
### Executive lines 182-215
ListResults<-function(fit.list,fit.sens,data.list,data.sens){
coef.list<-coef(fit.list )
vcov.list<-vcov( fit.list)

coef.sens<-coef(fit.sens)
vcov.sens<-vcov( fit.sens)

X.list<-as.matrix(cbind(1,na.omit(data.list[,3:ncol(data.list)])))
X.sens<-as.matrix(cbind(1,na.omit(data.list[,3:ncol(data.sens)])))

k<-length(coef(fit.list))/2
N<-1000

set.seed(11232)
coef.list.sim<-mvrnorm(N,coef.list[1:k],vcov.list[1:k,1:k])

set.seed(11232)
coef.sens.sim<-mvrnorm(N,coef.sens,vcov.sens)

#######
logistic <- function(x){exp(x)/(1+exp(x))}

results.list<-results.sens<-results.diff<-rep(NA,N)
for (i in 1:N){
    results.list[i]<-mean(logistic(X.list%*%coef.list.sim[i,]))
    results.sens[i]<-mean(logistic(X.sens%*%coef.sens.sim[i,]))
    results.diff[i]<-results.list[i]-results.sens[i]
  }

rbind(c(mean(results.list),quantile(results.list,c(0.025,0.975))),
        c(mean(results.sens),quantile(results.sens,c(0.025,0.975))),
        c(mean(results.diff),quantile(results.diff,c(0.025,0.975)))
  )
}
#########################################################################
# Test for Design Effect
print(de.test1)
print(de.test2)
print(de.test3)
print(de.test4)

##################### Results
results.list.1<-ListResults(fit.list.1,fit.sens.1,data.list.1,data.sens.1)
results.list.2<-ListResults(fit.list.2,fit.sens.2,data.list.2,data.sens.2)
results.list.3<-ListResults(fit.list.3,fit.sens.3,data.list.3,data.sens.3)
results.list.4<-ListResults(fit.list.4,fit.sens.4,data.list.4,data.sens.4)

cnames<-c("List","Direct","Diff")
rownames(results.list.1)<-cnames
rownames(results.list.2)<-cnames
rownames(results.list.3)<-cnames
rownames(results.list.4)<-cnames

print(results.list.1)
print(results.list.2)
print(results.list.3)
print(results.list.4)

################################ Figure 1 ########################################
#excute lines 243-274 
png(file="ListExperiments.png",
    width=3.5,height=2.5,units="in",res=900)
par(mfcol = c(1, 1), mar=c(1.35,5,0.5,0.1), cex = .55, cex.axis = .8)

plot(c(0.1,0.2,0.3), results.list.1[,1], type="n",
     pch = c(16,17,15), xlim = c(0, 1.9), ylim = c(-0.05,0.8), 
     xlab = "", ylab = "Estimated Proportion / Difference in Proportions",
     main = "", axes =F, cex.lab=.8)

abline(h = 0, lty = 2,lwd=0.8)
abline(h=results.list.1[3,1],lty=2,col="gray",lwd=0.8)

points(c(0.1, 0.2, 0.3), results.list.1[,1], pch = c(16,17,15))
points(c(0.6, 0.7, 0.8), results.list.2[,1], pch = c(16,17,15))
points(c(1.1, 1.2, 1.3), results.list.3[,1], pch = c(16,17,15))
points(c(1.6, 1.7, 1.8), results.list.4[,1], pch = c(16,17,15))

axis(side = 2, at = seq(from = 0, to = 1, by = 0.2))
axis(side = 1, at = c(0.2, 0.7, 1.2, 1.7), tick = FALSE,
     labels = c("Distrust Central Leaders", "Openly Criticizing Govt", 
                "Witnessed Corruption","Bribed Cadres"), cex.axis = 0.65, mgp = c(0,0,0))

segments(c(0.1, 0.2, 0.3),results.list.1[,2],c(0.1, 0.2, 0.3),results.list.1[,3],lwd=0.8)
segments(c(0.6, 0.7, 0.8),results.list.2[,2],c(0.6, 0.7, 0.8),results.list.2[,3],lwd=0.8)
segments(c(1.1, 1.2, 1.3),results.list.3[,2],c(1.1, 1.2, 1.3),results.list.3[,3],lwd=0.8)
segments(c(1.6, 1.7, 1.8),results.list.4[,2],c(1.6, 1.7, 1.8),results.list.4[,3],lwd=0.8)


legend(0.1,.8,c("List","Direct","Diff = List - Direct"),
       lty=1,pch=c(16,17,15),cex=.8,lwd=0.6,
       box.lwd=1,pt.cex=0.6,box.col = NA, bty="n")
dev.off()

